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ABSTRACT 



s 



We review mathematical methods for the treatment of a system of Bose particles with nonuniform density. 
Common pitfalls are pointed out. 



> 

o . 

^ ! 1 Introduction 

O, 

\^ ' The observation of Bose-Einstein condensation in atomic traps [l]-[3] has renewed theoretical interest in 

ON ■ studying a system of Bose particles with a nonuniform density. The experiments involve atoms of mass m 

condensed in a harmonic trap characterized by frequency u>. The relevant parameters, as discussed more 
fully in [4] and [5], have the following orders of magnitude: 



. a ~ 10 6 cm, 

£ I A = (2irh 2 /mkT) 1/2 ~ 10" 4 cm, 

L = (h/muj) 1/2 ~ 10~ 4 cm, 

na 3 ~ 10~ 3 , (1) 

where a is the S-wave scattering length, A the thermal wavelength, L the size of the system, and n is 
the average particle density. Thus, na 3 <C 1, a/L <C 1, a/A <C 1, and we can treat the system as a 
low-temperature weakly- interacting dilute gas using well-known methods [6]. 

In this paper, we review the following theoretical methods, and point out some common errors and 
misconceptions. 

(a) The pseudopotential method: One replaces the potential by a <5-function, and, with the help of a 
Bogolubov transformation, obtains a weak-coupling expansion. This method seems adequate to describe 
current experiments. 

(b) The self-consistent field method: This is variational, and capable of treating strongly-interacting 
cases in principle. It reduces to the pseudopotential method in appropriate limits; but one should not use 
a ^-function potential in a variational calculation, because such a potential in 3 spatial dimension has no 
effect, if treated exactly. 

(c) The Gaussian variational method: One uses a Gaussian trial wave function, and obtains results 
similar to the self-consistent field method. It has the advantage of easily adaptable to the description of 
time-dependent phenomena. 

All three methods are closely related to one another. Although we will not use Feynman graphs, it 
might be relevant to note that the Bogolubov transformation corresponds to summing "one-loop" graphs. 

1 Present address: Institute for Theoretical Atomic and Molecular Physics, Harvard-Smithsonian Center for Astrophysics, 
Cambridge, MA 02138. 
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The self-consistent field method consists of readjusting parameters in the one- loop sum in a variational 
sense. The Gaussian method includes all one-loop contribution, plus parts of higher-loop contributions. Of 
these methods, the straightforward one- loop summation yields a expansion in na 3 . On the other hand, the 
accuracy of the other methods cannot be ascertained. 

We do not attempt a comprehensive review, and the references cited are by no means complete. In a 
subject with such a long history, it is difficult to be aware of all relevant sources, and we apologize for any 
omissions [7]. 



2 Pseudopotential method 

In low-temperature calculations, one can replace the interatomic potential with a (^-function potential 
(47ra?i 2 /m)<5(r), where a > is the S-wave scattering length. Such a potential should be treated in low-order 
perturbation theory only, because the 5-function potential in 3 spatial dimensions has no effect if treated 
exactly [8]. 

Consider the Schrodinger equation 

-V 2 VXr)+. g( 5 3 (r)=^(r), (2) 

where E is 2ra/7i 2 times the energy. Taking the Fourier transform of both sides leads to 

(E-k 2 )0(k) =#(0), (3) 

where <f>(k) is the Fourier transform of tp(r). This must hold for all k, including k 2 = E. Therefore ip(0) = 0. 
This condition is automatically satisfied for all partial waves except S-waves. For S-waves, this seems to 
require that i/j(r) be discontinuous at r = 0; but because of the singular nature of the potential at r = 0, one 
should be careful. By treating the (5-function as the limit of a square well, one can verify that in 3 dimensions 
the S-waves are the same as those of a free particle except at r = 0, where they drop to zero discontinuously. 
This behavior, however, does not lead to any scattering, nor level shift. 

As a more careful analysis [9] shows, the (5-function potential must be regarded as the first term in a 
pseudopotential, which correctly reproduces all scattering phase shifts, and depends on these phase shifts as 
input parameters. For the simple case of a hard-sphere potential, for which there is only one parameter a, 
the hard-sphere diameter, the pseudopotential has the form 



^pseudo(r) 



Airaft 2 



S 3 (r) | r _^3 (r)v ^ r + 



(4) 



The differential operator (d/dr)r removes spurious singularities in the wave function, which would otherwise 
make the ground-state energy diverge. In low-order perturbation theory, one can replace this operator by a 
simple subtraction procedure for the energy. 

Using the first term in the pseudopotential, and making a Bogolubov transformation, one can calculate 
exactly the first few terms of an expansion for the ground-state energy per particle of a uniform hard-sphere 
Bose gas [10] [11]: 



En — 



2Trnah 2 



m 



128 fna 3 \ 1/2 „ /4tt 



15 V 7T J V 3 



Vsj na 3 \n(na 3 ) + 0(na 3 ) 



(5) 



This shows that the energy is not analytic at a = 0, and therefore cannot be continued to negative values of 
a. 

Fetter [12] has given a systematic treatment of the nonuniform case based on the Bogolubov transfor- 
mation. For application to the atomic-trap experiments, one can also adapt the results of the uniform case 
using a Thomas- Fermi approximation [5]. 

What happens when the scattering length is negative? In this case, the two-body scattering is dominated 
by the attractive part of the potential. The attraction may or may not be strong enough to produce a 
two-particle bound state; but it will lead to an iV-body bound state, for a sufficiently large N. The many- 
body problem is very different from that of the repulsive case, for it depends on the details of the potential. 
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Consider two potentials with the same negative scattering length, one everywhere attractive, and the other 
having a repulsive core. The many-body system with the purely attractive potential will collapse spatially, 
whereas the one with hard core will yield extensive energies. 

When the scattering length is negative, we must therefore use an actual potential, for example, one with 
a hard core plus an attractive tail. We can replace the hard-core part by a pseudopotential, and do the 
Bogolubov transformation. In this manner, it has been shown [6] [13] that a dilute Bose system with such a 
potential exhibits a first-order phase transition separating a gas phase from a liquid phase. The transition 
line terminates in a critical point. The Bose-Einstein condensation occurs in the liquid phase, and divides 
it into liquid I and liquid II. In short, one reproduces qualitatively the phase diagram of 4 He. Since the 
calculation is perturbative, one can follow the collapse of the phase diagram to that of the ideal Bose gas 
when the potential is turned off. 

The results described above pertains to an infinite Bose system in the thermodynamic limit, and may be 
modified for a finite system. The authors of Ref . [2] observed Bose condensation in an atomic trap filled with 
in 7 Li, which has a negative scattering length. In this case, the smallness of the system may have rendered 
the first-order transition indistinct. It is also possible that the observed phenomenon is metastable. 



3 Variational methods and the energy gap 

A variational approach to the ground-state energy of a Bose gas was introduced by Girardeau and Arnowitt 
[14], who use (wisely) an actual potential instead of the 5-function. Their results are for a uniform gas, 
but otherwise very similar to what we obtain later in the self-consistent field method. The excitation 
spectrum in this approach, however, contains an energy gap, contrary to general expectations [15], and to 
the perturbation-theory results in the hard-sphere case [10]. 

A variational method is designed to yield the best ground-state energy, and may not yield the excitation 
spectrum accurately. In this instance, we know that the gap should have been filled with phonons, which 
can be described using general principles. We may therefore view the variational principle as a way to obtain 
an approximate ground-state wave function, and build the phonon states upon this ground state. 

Feynman [16] argues that the low-lying excited states of a Bose system are purely phonon states, owing 
to the statistics, and suggests the following one-phonon wave function of momentum k: 

N 

* k (r 1 ,...,r w ) = ^e ikr ^ (ri,...,r w ), (6) 
where is the ground-state wave function. The excitation energy is 

where Sk is the liquid structure factor 

5 k = -(vf ia t eki*o>., (8) 

n 

where 

a = n- 1 /2^; a t +kap> (9) 
p 

with Ok the annihilation operator of a free particle of momentum Kk, and f2 the volume. From general 
principles [17] [18], we expect 

S^^, (10) 
k^o 2mc 

where c is the sound velocity as computed from the compressibility of the ground state. This leads to the 
phonon spectrum 

E k = hck. (11) 
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Feynman's argument can be sharpened by recasting it as the statement that the longitudinal sum-rules are 
saturated by the phonon states [17] [18]. All the equations above are verified in the dilute hard-sphere Bose 
gas [10] [18]. 

Other types of states that can be built upon * are those describing superfluid flow [10] [18]: 

N 

* SU pcrflow(ri, ,r N ) = J| e^^oOi, ...,r N ), (12) 

3 

where a(r) is the superfluid phase, which is related to the superfluid velocity v s through 

v s (r) = -Va(r). (13) 
m 

With this wave function one can describe quantized vortices. 



4 Self-consistent field method: ground state 

The self-consistent field method was first used by Bogolubov in his treatment of superconductivity, and is 
therefore also known as the "Bogolubov method." We shall follow a concise formulation due to De Gennes 
[19]. A recent review of similar methods is given by Griffin [20]. 

We use a quantized-field representation with field operator ^(x), and canonical conjugate i^^x): 

[*(x),*t(y)] =( j3 (x _ y)- (14) 

The Hamiltonian is 

H = j d 3 a;*(x) t /i(x)1'(x) + i J d 3 a;d 3 y*(y) t *(x) t y(x, y)tf(x)*(y), (15) 

where 

h 2 

fc(x) = - — V 2 + Kxt(x)-/i, (16) 
2m 

where V cxt (x) is the external potential, and fi is the chemical potential. We displace the field operator by 
its ground-state expectation: 

*(x)=V(x)+4>(x), (17) 
where V( x ) is the displaced field operator, and 

0(x) = (*(x)>. (18) 

The function </>(x) describes a nonuniform condensate, and ip(x) annihilates particles not in the condensate. 

Variational trial states are taken to be the eigenstates of an effective Hamiltonian that is quadratic in 
the field operators, chosen to be a mean-field version of the actual Hamiltonian. There are two variational 
functions p(x,y) and A(x, y), which will turn out to be p(x, y) = (x)ip(y)} , and A(x, y) = (^(x)^(y)). 
That is, p and A are, respectively, the direct and off-diagonal density correlation functions of the uncondensed 
particles. 

The effective Hamiltonian is chosen to be 

H cfi = + + ffWt + H (2 \ (19) 
where is independent of tp, and and terms are respectively linear and quadratic in ip'- 

- J d 3 xr(x)fe(x)0(x) + iy d 3 xd 3 #*(y)0*(x)T/(x,y)0(x)0(y), (20) 

#W = J d 3 xi/' t (x)/i(x)0(x) 

+ J d 3 xd 3 yV(x, y) V f (x) [0(y)p(y, x) + 0(x)p(y, y) + <f (y)I>(x, y)] , (21) 
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+ J d 3 xd 3 yV(x,y) [V t (y)V(x)^(x,y)+V t (x)V(x)ii(y,y)] 

+ l - J d 3 xd 3 yV(x,y) [V t (x)V t (y)r>(x,y) + V(x)^(y)£>*(x,y)] , (22) 

where 

i?(x,y) = p(x,y) + 0*(x)0(y), 

£>(x,y) ee A(x,y)+0(x)0(y). (23) 

For arbitrary p and A, we can diagonalize the effective Hamiltonian through a generalized Bogolubov 
transformation: 

V(x) = Y, [^K - • (24) 

n 

where r?„ and 77* are annihilation and creation operators satisfying 

[im'jln] — &nm- (25) 

To preserve the canonical commutators for ip(x), we require 

K(x)<(y) - <(xH(y)] = <5 3 (x - y). (26) 

n 

We now require H e g to have the form 

H cfi = Eq+Y CnVtiVn, (27) 
n 

or, equivalently, 

[H eS ,ri n ] = -e n T} n , 

[H eS ,nl] = e n r,l (28) 

These conditions determine the functions v n , u n and (j>. The condition for 4> comes from the requirement 
that in H e g the terms linear in vanish. The resulting equations are 

Mx)0(x) + J d 3 yV(x, y) [<f>(y)p(y, x) + 0(x)p(y, y) 

+ ^(y)A(x,y) + |0(y)| 2 0(x)] =0, (29) 

e n u n (x) = h(x)u„(x) + y <2 3 yF(x, y) [u„(x)i?(y, y) + u„(y)i?(y, x) - u n (y)D(x, y)] , (30) 

- e n v n (x) = fc(x)v n (x) + J d 3 yV(x,y)[v n (x)R(y,y) + v n (y)R(y,x) - u n (y)D'(x,y)] . (31) 

To determine p and A, let |$) be the lowest eigenstate of i? e ff : 

ff cff |$) =Eo\$). (32) 

We require that the ground-state energy be at a minimum with respect to variations in p and A: 

6(H) = 0, (33) 

where the brackets denote expectation value with respect to |$). In calculating the expectation value above, 
we use Wick's theorem: 

(^(^(yMxMy)) = (^ t (x)V(y))^ t (y)^(x)) + (^t(x) ? A(x))(^( y ) 7A ( y )) 

+<V> t (x)V> t (y))<V(x)V(y)>. (34) 



The calculation is facilitated by noting that S(H e s) = 0, since <£> is an eigenstate. Thus we can use the 
equivalent condition 

S(H - H eB ) = 0. (35) 

The results are as follows: 

p(x,y) = (^(xMy)) =$>„(x)<(y), 

71 

A(x,y) = (V(x)^(y))=-X;«n(xK(y). (36) 

n 

When these are substituted into (^9|). (pfl), and (|3l|), we have a system of coupled nonlinear integro-differential 
equations. 

5 The uniform case 

We put V cx t = 0, assume V(x, y) = V(x — y), and seek solutions with uniform density by putting 

<H r ) = 00, 
u k (r) = fi-^e^-'sfc, 

Wk(r) = n-^e^'i/k, (37) 

where SEfe, and </>o can be taken to be real. The condition ( p6|) becomes x\ — y\ = 1, which can be satisfied 
by putting 

x fe = cosh o-fc, 

t/ fc = sinh a k . (38) 

This parametrization simplifies the calculations. For example, the liquid-structure factor of the ground state 
can be represented in the form 

Sk = H — - (cosh 2tTfc — sinh 2a k — 1) 

n 

^ [sinh2<Tp sinh2o-| k+p | + sinh 2 a p sinh 2 C| p +k|] j (39) 
p 

where no is the condensate density given through 

n = n + r^ 1 sinh 2 cr p . (40) 
p 

The Bogolubov result is recovered by neglecting the p sums, which is equivalent to putting no ~ n: 

1/2 



s k = 



(41) 



_e k + 2nV k _ 

where e k = h 2 k 2 /2m, and V k is the Fourier transform of the interaction potential. 

6 Reduction to ^-function potential 

The equations in the last section are rather complex, and difficult to solve even numerically. To gain some 
insight, we reduce them to a simpler form by using a ^-function potential. We shall show that this potential 
has no effect if self-consistency is enforced, at least in the case of a uniform system. We put 

F(x,y) = ff <5 3 (x-y), (42) 



G 



where 

g = 4irh 2 a/m. (43) 

The functions p and A now depend only on x: 

p(x) = (^(xWx))=^„(x)<(x) 

/I 

A(x) - (V(x)V>(x))=-5>„(x)<(x). (44) 

n 

Note that p is the depletion of the condensate due to interactions. The equations (p9|), (|3C|), and ( pl| ) reduce 
to 

[/i + 2.gp + .g|0| 2 ] ( / ) + .gA0*=O, (45) 

(h + 2gp + 2gcf)*(j>)u n ~ g{A + (f) 2 )v n = eu„, 
(ft + 2.gp + 2gcj)*(j))v n - g(A* + <j)* 2 )u n = -ev n , (46) 

where we have suppressed the x dependence of the functions. The dilute uniform hard-sphere gas is recovered 
by setting p = A = 0, for they give higher-order contributions. 



For p — A — 0, (45) is the "nonlinear Schrodinger equation" [21] [22], which has been widely used in 
numerical studies of the condensate and its excitations [23]. One should be aware, however, that p and A 
may be important, especially for the excitations. 

In the uniform case, with V^xt = 0, we obtain, in the notation of the last section, 

tanh 2a k = ^ ' >— , (47) 

e k + g(n- p- A) 

where e k = fi 2 k 2 /2m, and the chemical potential p has been eliminated in terms of the density n. The 



self-consistency conditions (44) then lead to 



P = 4* 



l r 

^ / dkk 2 {f k [e k + g{n - p - A)] - 1} , 
A = dkk 2 f k (n-p + A), (48) 

where 

h = [e 2 k + 2ge k (n - p — A) — 4g 2 (n - p)A] ~ 1/2 . (49) 

The integrals are cut off at A, because otherwise A would be divergent. Solving for A from the second 
equation, we obtain in the limit A — > oo 

A - p + n = 0. (50) 
The first equation then gives, in the limit A — > oo, 

p = Q. (51) 

That is, there is no depletion of the condensate. It is then straightforward to show that the system is an ideal 
Bose gas. This result is reassuring; it shows that the method is able to verify that the 5-function potential 
has no effect. 

7 Self-consistent field method: finite temperatures 

To extend the variational approach to finite temperatures, we can use the eigenstates of H g to calculate 
the partition function, and then minimize the free energy 

F=(H) - TS, (52) 
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where the brackets now denote thermal average with respect to H c g, and S is the entropy. De Gennes [19] 
argues that the quantity 

F cS = (H cS )-TS (53) 

is stationary, and implements the variational principle by minimizing F— F e s . But the assertion is correct only 
if S were defined with respect to the eigenvalue spectrum of H c g . In fact, as we explain below, S is defined 
with respect to the energy levels (a\H\a), where \a) are the eigenstates of H c g. The finite-temperature 
results in Ref.[19] are therefore incorrect. 

A proper variational principle for the free energy [24] is based on the inequality 

£<a|e-H«>>£e-'WH (54) 

a a 

This shows that, to calculate the trial partition function, we must use the energy spectrum 

W a = (a\H\a). (55) 

In contrast, the other scheme leads to a free Bose gas with single-particle energy e„, the eigenvalues from 
( |30| ) and ( |3lf ). Because of the simplicity, it is widely used in the literature; but it does not follow from a 
variational principle. 

As an example of the difference between the two schemes, the calculation of the partition of a hard-sphere 
Bose gas in Ref.[25] corresponds to using (|55|), and yields the virial coefficients in the gas phase. On the other 
hand, a calculation based on the other scheme would give a free Bose gas above the transition temperature. 

8 Gaussian variational method 

We go to a representation in which the field operator ^(x) and its conjugate have the forms 

1 



•w - ^ 

»'<«) = £ 



yj(x) 



y>(x) - 



6 

6<p(x) 



(56) 



where y(x) is a c-number function. The more familiar representation, in which '5 is diagonal, would be 
awkward to use in a non-relativistic case, in which \& and are canoncial conjugates. (A similar situation 
happens in a relativistic fermion field obeying the Dirac equation.) 

The state of the system is described by a wave functional $[<£>]. We use a Gaussian form for a trial wave 
functional: 

<%]=Cexpj-^ J d 3 xd 3 2 /[^(x)-^o(x)]G- 1 (x,y)Ky)-^ (y)]|, (57) 

where G is a normalization constant and ipo and G(x, y) are the variational parameters. Expectation values 
are given by functional integrals, for example 



(H)= J (D<p)&*[<p]H$[<p] (58) 

where $ is normalized to unity. It is easy to show that 

(*(x)> = ^o(x), 
MxMy)) = G(x,y)+^ (x)^o(y). (59) 



The variational principle states that 



<5G(x,y) 6<p (x) 
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The equation obtained for G and ipo will not be given directly. To make contact with the self-consistent 
field method, we quote the results 



p(x,y) ee (V>(x)V(y)) = 



A(x,y) 



Mx)V(y)) 



f 

2 

1 

~2 



G- 1 (x,y)+G(x,y)- ( 5 3 (x-y) 



G- 1 (x,y)~G(x,y) 



(61) 



For further analysis, it is convenient to introduce 

1 



*(x,y) = 



2V2 
1 

5\/2 



G- 1 / 2 (x,y)+2G 1 /2 (x , y ) 
G- 1 / 2 (x,y)-2G 1 /2 (x , y ) 



(62) 



it easy to show that p(x, y) = J d 3 zF(x, z)Y (z, y), and A(x, y) = — J d 3 zX(x,z)Y(z,y). Now expand X 
and Y in terms of a basis {6„(x)}: 



X(x,y)=^X„MxR(y), 

n 

y(x,y)=5>„&„(x)&:(y), 

n 

where X n and Y n are real. We then find 

P (x,y) = ^6„(x)6;(y)F„ 2 , 

A(x,y) = -£& n (x)&*(y)A* 

n 

Comparison with ( |36| ) leads to the identification 

U„(x) = 6 n (x)X„, 
u n (x) = 6„(x)F„. 



(63) 



(64) 



(65) 



This shows the equivalence between this method and the self-consistent field method. 

One of the advantages of the Gaussian method lies in the possibility of generalization to the "time- 
dependent Hartree-Fock" method [26] [27], which has been used successfully in nuclear physics. One of us 
(P.T.) plans to discuss this in a separate publication. 
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